﻿using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Windows.Forms;

namespace Lab
{
    public partial class Form1 : Form
    {
        public Form1()
        {
            InitializeComponent();
        }
        int n = 2; // Число разбиений по X
        int m = 2; // Число разбиений по Y
        double h = 0.5;
        double k = 0.5;
        double[] NetX;
        double[] NetY;
        double w;
        double[,] L;
        double[,] D;
        double[,] R;
        double[,] A;
        double[] b;
        double[] v_next;
        int size;
        int num_steps = 0;
        double Epsil;

        double U(double x, double y)
        {
            return Math.Sin(Math.PI * x * y);
        }

        double f(double x, double y)
        {
            return -Math.Sin(Math.PI * x * y) * (Math.Pow(Math.PI * x, 2.0) + Math.Pow(Math.PI * y, 2.0));
        }

        double[,] SetA()
        {
            double[,] res = new double[size, size];
            int count = 0;
            //Заполним диагональ и правую матрицу
            for (int i = 1; i < size; i++)
            {
                res[i - 1,i - 1] = -2.0 * (1.0 / Math.Pow(h, 2.0) + 1.0 / Math.Pow(k, 2.0));
                if (count != n - 2)
                {
                    res[i - 1, i] = 1.0 / Math.Pow(h, 2.0);
                    count++;
                }
                else
                {
                    res[i - 1, i] = 0.0;
                    count = 0;
                }
                if (i >= n - 1)
                {
                    res[i - n + 1, i] = 1.0 / Math.Pow(k, 2.0);
                }
                if(i == size - 1)
                    res[i, i] = -2.0 * (1.0 / Math.Pow(h, 2.0) + 1.0 / Math.Pow(k, 2.0));
            }
            //Заполним левую матрицу
            count = 0;
            for (int j = 1; j < size; j++)
            {
                if (count != n - 2)
                {
                    res[j, j - 1] = 1.0 / Math.Pow(h, 2.0);
                    count++;
                }
                else
                {
                    res[j, j - 1] = 0.0;
                    count = 0;
                }
                if (j >= n - 1)
                {
                    res[j, j - n + 1] = 1.0 / Math.Pow(k, 2.0);
                }
            }
            return res;
        }        

        double[,] SetL()
        {
            double[,] res = new double[size, size];
            int count = 0;
            for (int j = 1; j < size; j++)
            {
                if (count != n - 2)
                {
                    res[j, j - 1] = 1.0 / Math.Pow(h, 2.0);
                    count++;
                }
                else
                {
                    res[j, j - 1] = 0.0;
                    count = 0;
                }
                if (j >= n - 1)
                {
                    res[j, j - n + 1] = 1.0 / Math.Pow(k, 2.0);
                }
            }
            return res;
        }

        double[,] SetD()
        {
            double[,] res = new double[size, size];
            for (int i = 0; i < size; i++)
            {
                res[i, i] = -2.0 * (1.0 / Math.Pow(h, 2.0) + 1.0 / Math.Pow(k, 2.0));
            }
            return res;
        }

        double[,] SetR()
        {
            double[,] res = new double[size, size];
            int count = 0;            
            for (int i = 1; i < size; i++)
            {                
                if (count != n - 2)
                {
                    res[i - 1, i] = 1.0 / Math.Pow(h, 2.0);
                    count++;
                }
                else
                {
                    res[i - 1, i] = 0.0;
                    count = 0;
                }
                if (i >= n - 1)
                {
                    res[i - n + 1, i] = 1.0 / Math.Pow(k, 2.0);
                }             
            }
            return res;
        }

        double[] Setb()
        {
            double[] res = new double[size];
            int t = 1;
            for (int i = 0; i < size; i++)
            {
                if (i < n - 1)
                {
                    if (i == 0)
                        res[i] = f(NetX[i + 1], NetY[t]) - (1.0 / Math.Pow(h, 2.0)) * U(NetX[i], NetY[t]) - (1.0 / Math.Pow(k, 2.0)) * U(NetX[i + 1], NetY[i]);
                    else
                    {
                        if(i == n - 2)
                            res[i] = f(NetX[i + 1], NetY[t]) - (1.0 / Math.Pow(h, 2.0)) * U(NetX[n], NetY[t]) - (1.0 / Math.Pow(k, 2.0)) * U(NetX[i + 1], NetY[0]);
                        else
                        {
                            res[i] = f(NetX[i + 1], NetY[t]) - (1.0 / Math.Pow(k, 2.0)) * U(NetX[i + 1], NetY[0]);
                        }
                    }
                }
                if(i % (n - 1) == 0 && i != 0)
                    t++;
                if (i >= n - 1 && i < size - n + 1)
                {
                    if (i % (n - 1) == 0)
                    {
                        res[i] = f(NetX[i - (t - 1) * (n - 1) + 1], NetY[t]) - (1.0 / Math.Pow(h, 2.0)) * U(NetX[0], NetY[t]);
                    }
                    else
                    {
                        if ((i - 2) % (n - 1) == 0)
                        {
                            res[i] = f(NetX[i - (t - 1) * (n - 1) + 1], NetY[t]) - (1.0 / Math.Pow(h, 2.0)) * U(NetX[n], NetY[t]);
                        }
                        else
                        {
                            res[i] = f(NetX[i - (t - 1) * (n - 1) + 1], NetY[t]);
                        }
                    }
                }
                if (i >= size - n + 1)
                {
                    if (i == size - n + 1)
                    {
                        res[i] = f(NetX[i - (t - 1) * (n - 1) + 1], NetY[t]) - (1.0 / Math.Pow(h, 2.0)) * U(NetX[0], NetY[t]) - (1.0 / Math.Pow(k, 2.0)) * U(NetX[i - (t - 1) * (n - 1) + 1], NetY[m]);
                    }
                    else
                    {
                        if (i == size - 1)
                        {
                            res[i] = f(NetX[i - (t - 1) * (n - 1) + 1], NetY[t]) - (1.0 / Math.Pow(h, 2.0)) * U(NetX[n], NetY[t]) - (1.0 / Math.Pow(k, 2.0)) * U(NetX[i - (t - 1) * (n - 1) + 1], NetY[m]);
                        }
                        else
                        {
                            res[i] = f(NetX[i - (t - 1) * (n - 1) + 1], NetY[t]) - (1.0 / Math.Pow(k, 2.0)) * U(NetX[i - (t - 1) * (n - 1) + 1], NetY[m]);
                        }
                    }
                }
            }
            return res;
        }

        double SetW()
        {
            double res = 0.0;
            double max = 0.0;
            if (n != 2 && m != 2)
            {
                double[,] B = new double[size, size];
                double[,] Diag = new double[size, size];
                Diag = D;
                double[,] LR = new double[size, size];
                for (int i = 0; i < size; i++)
                {
                    for (int j = 0; j < size; j++)
                    {
                        LR[i, j] = L[i, j] + R[i, j];
                    }
                }
                int info;
                alglib.matinvreport rep;
                alglib.rmatrixinverse(ref Diag, out info, out rep);
                double tmp = 0.0;
                for (int i = 0; i < size; i++)
                {
                    for (int j = 0; j < size; j++)
                    {
                        for (int z = 0; z < size; z++)
                        {
                            tmp = tmp + Diag[i, z] * LR[z, j];
                        }
                        B[i, j] = tmp;
                        tmp = 0.0;
                    }
                }
                double[] evd = new double[size];
                double[] q;
                double[,] t, r;
                alglib.rmatrixevd(B, size, 0, out evd, out q, out t, out r);
                max = Math.Abs(evd[0]);
                for (int i = 1; i < size; i++)
                {
                    if (max < Math.Abs(evd[i]))
                        max = Math.Abs(evd[i]);
                }
                /*for (int i = 0; i < size; i++)
                {
                    max = max + Math.Abs(B[0, i]);
                }
                max = max - Math.Abs(B[0, 0]);
                double tmp_max = 0.0;
                for (int i = 1; i < size; i++)
                {
                    for (int j = 0; j < size; j++)
                    {
                        tmp_max = tmp_max + Math.Abs(B[i, j]);
                    }
                    tmp_max = tmp_max - Math.Abs(B[i, i]);
                    if (max < tmp_max)
                    {
                        max = tmp_max;
                        tmp_max = 0.0;
                    }
                    else
                    {
                        tmp_max = 0.0;
                    }
                }*/
            }
            else
            {
                max = 0.0;
            }
            res = 2.0 / (1.0 + Math.Sqrt(1.0 - Math.Pow(max, 2.0)));
            return res;
        }

        void Spravka(double epsilon_max, double tochnost)
        {
            label9.Text = "(m, n) = (" + Convert.ToString(n) + ", " + Convert.ToString(m) + ").";
            label11.Text = "Точность метода: " + Eps.Text;
            label12.Text = "Параметр метода: " + Convert.ToString(w);
            label13.Text = "Количество итераций на решение: " + Convert.ToString(num_steps);
            label14.Text = "Достигнутая точность метода: " + Convert.ToString(epsilon_max);
            label15.Text = "Задача решена с точностью: " + Convert.ToString(tochnost);
        }

        private void btnSolve_Click(object sender, EventArgs e)
        {
            n = Convert.ToInt32(numX.Value);
            m = Convert.ToInt32(numY.Value);
            size = n * m - n - m + 1;
            h = 1.0 / n;
            k = 1.0 / m;
            NetX = new double[n + 1];
            NetY = new double[m + 1];
            NetX[0] = 1.0;
            for (int i = 1; i < n + 1; i++)
            {
                NetX[i] = NetX[i - 1] + h;
            }
            NetY[0] = 2.0;
            for (int j = 1; j < m + 1; j++)
            {
                NetY[j] = NetY[j - 1] + k;
            }            
            v_next = new double[size];
            Epsil = Convert.ToDouble(Eps.Text);
            A = new double[size, size];
            L = new double[size, size];
            D = new double[size, size];
            R = new double[size, size];
            b = new double[size];
            A = SetA();
            L = SetL();
            D = SetD();
            R = SetR();
            b = Setb();
            w = SetW();
            double eps_max = 0.0; // текущее значение прироста
            double MaxToch = 0.0; //Точность
            if (n == 2 && m == 2)
            {
                double cur_v = 0.0;
                v_next[0] = (1.0 / (-2.0 / Math.Pow(h, 2.0) - 2.0 / Math.Pow(k, 2.0))) * (f(NetX[1], NetY[1]) - U(NetX[0], NetY[1]) / Math.Pow(h, 2.0) - U(NetX[2], NetY[1]) / Math.Pow(h, 2.0) - U(NetX[1], NetY[0]) / Math.Pow(k, 2.0) - U(NetX[1], NetY[2]) / Math.Pow(k, 2.0));
                while (Math.Abs(cur_v - v_next[0]) > Epsil)
                {
                    cur_v = v_next[0];
                    v_next[0] = (1.0 / (-2.0 / Math.Pow(h, 2.0) - 2.0 / Math.Pow(k, 2.0))) * ((1 - w) * A[0,0] * cur_v + w * b[0]);
                }
                eps_max = Math.Abs(cur_v - v_next[0]);
                MaxToch = Math.Abs(U(NetX[1], NetY[1]) - v_next[0]);
            }
            else
            {                
                double eps_cur = 0.0; // для подсчета текущего значения прироста
                double x_old; // старое значение преобразуемой компоненты вектора v_next
                double x_new; // новое значение преобразуемой компоненты вектора v_next
                bool flag = false; //условие остановки
                while (!flag)
                {
                    eps_max = 0.0;
                    for (int i = 0; i < size; i++)
                    {
                        x_old = v_next[i];
                        x_new = (1.0 - w) * A[i, i] * v_next[i] + w * b[i];
                        for (int j = 0; j < size; j++)
                        {
                            if (j != i)
                            {
                                x_new = x_new - w * A[i, j] * v_next[j];
                            }
                        }
                        x_new = x_new / A[i, i];
                        eps_cur = Math.Abs(x_old - x_new);
                        if (eps_cur > eps_max)
                        {
                            eps_max = eps_cur;
                        }
                        v_next[i] = x_new;
                    }
                    num_steps++;
                    if ((eps_max < Epsil) || (num_steps >= Convert.ToInt32(numSTEPS.Value)))
                        flag = true;
                }
                MaxToch = Math.Abs(U(NetX[1], NetY[1]) - v_next[0]);
                for (int j = 1; j < m; j++)
                {
                    for (int i = 1; i < n; i++)
                    {
                        if (MaxToch <= Math.Abs(U(NetX[i], NetY[j]) - v_next[i - 1 + (j - 1) * (n - 1)]))
                        {
                            MaxToch = Math.Abs(U(NetX[i], NetY[j]) - v_next[i - 1 + (j - 1) * (n - 1)]);
                        }
                    }
                }
            }

            if (dataGridView1.Rows.Count != 1)
            {
                dataGridView1.Rows.Clear();
            }
            if (dataGridView1.Columns.Count != 1)
            {
                dataGridView1.Columns.Clear();
            }
            dataGridView1.Columns.Add("i,j", "i,j");
            for (int i = 0; i < n + 1; i++)
            {
                dataGridView1.Columns.Add(Convert.ToString(i), "x[" + Convert.ToString(i) + "] = " + Convert.ToString(NetX[i]));
            }
            dataGridView1.Rows.Add(m + 1);
            for (int j = 0; j < m + 1; j++)
            {
                dataGridView1.Rows[j].Cells[0].Value = "y[" + Convert.ToString(j) + "] = " + Convert.ToString(NetY[j]);
                dataGridView1.Rows[j].Cells[1].Value = Convert.ToString(U(NetX[0], NetY[j]));
            }
            for (int i = 2; i < n + 2; i++)
            {
                dataGridView1.Rows[0].Cells[i].Value = Convert.ToString(U(NetX[i - 1], NetY[0]));
            }
            for (int i = 2; i < n + 2; i++)
            {
                dataGridView1.Rows[m].Cells[i].Value = Convert.ToString(U(NetX[i - 1], NetY[m]));
            }
            for (int j = 1; j < m; j++)
            {
                dataGridView1.Rows[j].Cells[n + 1].Value = Convert.ToString(U(NetX[n], NetY[j]));
            }
            for (int j = 1; j < m; j++)
            {
                for (int i = 2; i < n + 1; i++)
                {
                    dataGridView1.Rows[j].Cells[i].Value = Convert.ToString(v_next[i - 2 + (j - 1) * (n - 1)]);
                }
            }        
            for (int j = 0; j < m + 1; j++)
            {
                dataGridView1.Rows[j].Cells[0].Value = "y[" + Convert.ToString(j) + "] = " + Convert.ToString(NetY[j]);
            }

            if (dataGridView2.Rows.Count != 1)
            {
                dataGridView2.Rows.Clear();
            }
            if (dataGridView2.Columns.Count != 1)
            {
                dataGridView2.Columns.Clear();
            }
            dataGridView2.Columns.Add("(i,j)", "(x, y)");
            for (int i = 0; i < n + 1; i++)
            {
                dataGridView2.Columns.Add(Convert.ToString(i), "x[" + Convert.ToString(i) + "] = " + Convert.ToString(NetX[i]));
            }
            dataGridView2.Rows.Add(m);
            for (int j = 0; j < m + 1; j++)
            {
                dataGridView2.Rows[j].Cells[0].Value = "y[" + Convert.ToString(j) + "] = " + Convert.ToString(NetY[j]);
            }
            for (int j = 0; j < m + 1; j++)
            {
                for (int i = 0; i < n + 1; i++)
                {
                    dataGridView2.Rows[j].Cells[i + 1].Value = Convert.ToString(U(NetX[i], NetY[j]));
                }
            }

            if (dataGridView3.Rows.Count != 1)
            {
                dataGridView3.Rows.Clear();
            }
            if (dataGridView3.Columns.Count != 1)
            {
                dataGridView3.Columns.Clear();
            }
            dataGridView3.Columns.Add("(x, y)", "(x, y)");
            for(int i = 0; i < n + 1; i++)
            {
                dataGridView3.Columns.Add(Convert.ToString(i), "x[" + Convert.ToString(i) + "] = " + Convert.ToString(NetX[i]));
            }
            dataGridView3.Rows.Add(m);
            for (int j = 0; j < m + 1; j++)
            {
                dataGridView3.Rows[j].Cells[0].Value = "y[" + Convert.ToString(j) + "] = " + Convert.ToString(NetY[j]);
            }
            for(int i = 0; i < n + 1; i++)
            {
                dataGridView3.Rows[0].Cells[i + 1].Value = "0";
            }
            for(int j = 0; j < m + 1; j++)
            {
                dataGridView3.Rows[j].Cells[1].Value = "0";
            }
            for(int i = 1; i < n + 1; i++)
            {
                dataGridView3.Rows[m].Cells[i].Value = "0";
            }
            for(int j = 1; j < m + 1; j++)
            {
                dataGridView3.Rows[j].Cells[n + 1].Value = "0";
            }
            for (int i = 2; i < n + 1; i++) 
            {
                for (int j = 1; j < m; j++)
                {
                    dataGridView3.Rows[j].Cells[i].Value = Convert.ToString(Math.Abs(U(NetX[i - 1], NetY[j]) - v_next[i - 2 + (j - 1) * (n - 1)]));
                }
            }
            Spravka(eps_max, MaxToch);
        }//btnSolve_Click

    }
}
